######################### Load Libraries #######################################

library(foreign)

library(ggplot2)

library(lattice)

library(plyr)

library(gridExtra)

############################ Load Data #########################################

# set working directory
setwd("~/_data/")

# load LPR data set
load("lprdata.RData")

################################################################################
########### Austria: LPR calculation using Silverman's Rule of Thumb ###########
################################################################################

###########################################
### calculate LPR using Gaussian kernel ###
###########################################

data$lpr_ga_nrd0 <- rep(NA,nrow(data))

for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="gaussian")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_ga_nrd0"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

###############################################
### calculate LPR using Epanechnikov kernel ###
###############################################

data$lpr_ep_nrd0 <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="epanechnikov")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_ep_nrd0"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

##############################################
### calculate LPR using rectangular kernel ###
##############################################

data$lpr_re_nrd0 <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="rectangular")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_re_nrd0"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

#############################################
### calculate LPR using triangular kernel ###
#############################################

data$lpr_tr_nrd0 <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="triangular")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_tr_nrd0"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

###########################################
### calculate LPR using biweight kernel ###
###########################################

data$lpr_bi_nrd0 <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="biweight")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_bi_nrd0"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

############################################
### calculate LPR using optcosine kernel ###
############################################

data$lpr_op_nrd0 <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="optcosine")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_op_nrd0"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

###########################################
### combine results into new data frame ###
###########################################

data <- data[data$isocode=="AUT",]

lpr_nrd0 <- matrix(NA,120,2)

lpr_nrd0[,1] <- as.numeric(c(data$lpr_ga_nrd0,data$lpr_ep_nrd0,data$lpr_re_nrd0,data$lpr_tr_nrd0,data$lpr_bi_nrd0,data$lpr_op_nrd0))

lpr_nrd0[,2] <- as.factor(c(rep("Gaussian",20),rep("Epanechnikov",20),rep("Rectangular",20),rep("Triangular",20),rep("Biweight",20),rep("Cosine",20)))

lpr_nrd0 <- as.data.frame(lpr_nrd0)

colnames(lpr_nrd0) <- c("LPR","Kernel")


#############################################################################
########### Austria: LPR calculation using Sheather & Jones' Pilot ##########
########### Estimation of Derivatives                              ##########
#############################################################################

###########################################
### calculate LPR using Gaussian kernel ###
###########################################

data$lpr_ga_sj <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="gaussian", bw = "SJ")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_ga_sj"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

###############################################
### calculate LPR using Epanechnikov kernel ###
###############################################

data$lpr_ep_sj <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="epanechnikov", bw = "SJ")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_ep_sj"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

##############################################
### calculate LPR using rectangular kernel ###
##############################################

data$lpr_re_sj <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="rectangular", bw = "SJ")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_re_sj"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

#############################################
### calculate LPR using triangular kernel ###
#############################################

data$lpr_tr_sj <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="triangular", bw = "SJ")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_tr_sj"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

###########################################
### calculate LPR using biweight kernel ###
###########################################

data$lpr_bi_sj <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="biweight", bw = "SJ")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_bi_sj"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

############################################
### calculate LPR using optcosine kernel ###
############################################

data$lpr_op_sj <- rep(NA,nrow(data))

# calculate LPR
for (i in 7:nrow(data[data$isocode=="AUT",])) {
	
	z.pr <- density(data[data$isocode=="AUT","s.pr"][-c(i:nrow(data[data$isocode=="AUT",]))], kernel="optcosine", bw = "SJ")
	f.pr <- approxfun(z.pr$x, z.pr$y, yleft = 0, yright = 0)
	data[data$isocode=="AUT","lpr_op_sj"][i] <- integrate(f.pr, -Inf, -data[data$isocode=="AUT","seatsharegap"][i])$value

}

###########################################
### combine results into new data frame ###
###########################################

data <- data[data$isocode=="AUT",]

lpr_sj <- matrix(NA,120,2)

lpr_sj[,1] <- as.numeric(c(data$lpr_ga_sj,data$lpr_ep_sj,data$lpr_re_sj,data$lpr_tr_sj,data$lpr_bi_sj,data$lpr_op_sj))

lpr_sj[,2] <- as.factor(c(rep("Gaussian",20),rep("Epanechnikov",20),rep("Rectangular",20),rep("Triangular",20),rep("Biweight",20),rep("Cosine",20)))

lpr_sj <- as.data.frame(lpr_sj)

colnames(lpr_sj) <- c("LPR","Kernel")

################################## SM Figure 4 #################################

# set working directory for SM graphs
setwd("~/_supportingmaterials_tabfig/")

# generate pdf graph (black and white)
pdf("aus_density_nrd.pdf")
bwplot(LPR ~ as.factor(Kernel), ylab="LPR (Austria)", xlab="Kernel",main="Silverman's Rule of Thumb",lpr_nrd0,par.settings = list(plot.symbol = list(col = "black"),box.umbrella = list(col ="black"), box.rectangle=list(col="black")),fill="light grey",scales=list(x=list(cex=0.8, labels=c("Gaussian","Epanechnikov","Rectangular","Triangular","Biweight","Cosine") )))
dev.off()

# generate pdf graph (black and white)
pdf("aus_density_sj.pdf")
bwplot(LPR ~ as.factor(Kernel), ylab="LPR (Austria)", xlab="Kernel",main="Sheather & Jones' Pilot Estimation of Derivatives",lpr_sj,par.settings = list(plot.symbol = list(col = "black"),box.umbrella = list(col ="black"), box.rectangle=list(col="black")),fill="light grey",scales=list(x=list(cex=0.8, labels=c("Gaussian","Epanechnikov","Rectangular","Triangular","Biweight","Cosine") )))
dev.off()